Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Arrayslice #96

Merged
merged 57 commits into from
Sep 29, 2017
Merged

Arrayslice #96

merged 57 commits into from
Sep 29, 2017

Conversation

HomerReid
Copy link
Contributor

@HomerReid HomerReid commented Sep 16, 2017

The primary new routine here is

 std::complex<double> *get_array_slice(const volume &where,
                                        std::vector<component> components,
                                        field_function fun=0, void *fun_data_=0,
                                        std::complex<double> *slice=0);

The implementation is similar to output_hdf5, but I rearranged the calling convention based on the following rationale:
-- the subvolume of interest, and the desired component(s), must always be specified and have no defaults, so those should come first
-- for components, I use a std::vector instead of pointer-plus-integer-valued-size-parameter as in output_hdf5
-- callers should be able to omit the field_function argument to request the identity function
-- callers should have the option of (a) passing a caller-allocated array to be filled in with the array-slice data, (b) having the function allocate and return a new array.

For option (a) above I added two new routines:
get_array_slice_dimensions()
and
allocate_array_slice_buffer()
with self-explanatory calling conventions.

If the caller passes a non-null buffer, it's not obvious how we would check to make sure its the correct size. Some options include
(a) don't check, creating the risk of core dump if the caller passes a buffer that's too small
(b) use a container like std::vector<cdouble> instead of a bare pointer
(c) use some sort of high-level multidimensional array construct that would play nicely with existing python libraries

For the time being I've gone with (a).

Testing: Added tests/array_slice_test.cpp for testing. This is a modified version of h5test.cpp. The implementation is not complete, but the idea is to use the same volumes and subvolumes as h5test.cpp and compare the resulting array slices with arrays read in from hdf5 files produced by running h5test.

Usage example:

 int dims[3], dirs[3];
 int rank=f.get_array_slice_dimensions(subvolume, dims, dirs);

 std::vector<component> components(Ez);
 components.push_back(Hz);
 cdouble *slice=f.get_array_slice(subvolume, components);

 int NC=components.size();
 int NX=dims[0];
 int NY=dims[1];
 for(int nx=0; nx<NX; nx++)
   for(int ny=0; ny<NY; ny++)
    printf("At (%i,%i): (Ez,Hz)=(%e,%e)\n",nx,ny,
               real(slice[0*NX*NY + nx*NY + ny]),
               real(slice[1*NX*NY + nx*NY + ny]));
```

@HomerReid
Copy link
Contributor Author

HomerReid commented Sep 26, 2017

The latest commit to this branch includes the following updates:

-- the array_slice suite of routines in src/array_slice.cpp, which includes routines like get_array_slice_dimensions and get_array_slice, is now working

-- added working unit test libmeepgeom/array-slice-ll

-- I also added a couple of typemaps in python/meep.i and added python/examples/cavity_arrayslice.py demonstrating manipulation of array slices as numpy arrays in python codes. This is the same geometry as the holey_wvg_cavity example. The python script obtains and plots 1D and 2D array slices after the fields have finished. Here is the relevant python snippet, from the bottom of cavity_arrayslice.py:

# 1D slice
dims1d=np.zeros(3,dtype=np.int32)
v1d=mp.volume( mp.vec(xMin, 0.0), mp.vec(xMax, 0.0) )
rank1d  = sim.fields.get_array_slice_dimensions(v1d, dims1d);
NX1 = dims1d[0];
slice1d = np.zeros(NX1, dtype=np.double);
sim.fields.get_array_slice(v1d, mp.Hz, slice1d);

# 2D slice
dims2d=np.zeros(3,dtype=np.int32)
v2d=mp.volume( mp.vec(xMin, yMin), mp.vec(xMax, yMax) )
rank2d  = sim.fields.get_array_slice_dimensions(v2d, dims2d);
NX2 = dims2d[0];
NY2 = dims2d[1];
N2 = NX2*NY2;
slice2d = np.zeros(N2, dtype=np.double);
sim.fields.get_array_slice(v2d, mp.Hz, slice2d);

# plot 1D slice
plt.subplot(1,2,1);
x1d=np.linspace(xMin, xMax, dims1d[0]);
plt.plot(x1d, slice1d);

# plot 2D slice
plt.subplot(1,2,2);
dy = (yMax-yMin) / dims2d[1];
dx = (xMax-xMin) / dims2d[0];
(x2d,y2d)=np.mgrid[ slice(xMin, xMax, dx), slice(yMin, yMax, dy)];
plt.contourf( x2d, y2d, np.reshape(slice2d, (dims2d[0], dims2d[1])) )
plt.colorbar();
plt.show()

Here's a screenshot of the output this produces:

screenshot

@oskooi @ChristopherHogan

@HomerReid
Copy link
Contributor Author

Just to update, all tests are passing on my local machines, and the individual tests all pass on travis, but the overall make check fails on travis due to the annoying re-linking error we saw before.

Previously it seems Chris fixed this by doing the python2 and python3 builds in different directories, but that doesn't seem to be working here.

Does anybody have any other ideas what might be going on?

I'll keep plugging away, but in the meantime the array_slice feature is ready to use per the python snippet above.

@oskooi @ChristopherHogan

@ChristopherHogan
Copy link
Contributor

Recently I've seen this error message when the real cause was a previously failed build step. I had to search through the logs for "error" and "fail." Looking at the travis log from your latest commit, I see this:

make[2]: *** No rule to make target `array_slice_test.o', needed by `array_slice_test'.  Stop.
make[2]: Leaving directory `/home/travis/build/stevengj/meep/build/tests'
make[1]: *** [check-am] Error 2
make[1]: Leaving directory `/home/travis/build/stevengj/meep/build/tests'
make: *** [check-recursive] Error 1

@HomerReid
Copy link
Contributor Author

Thanks Chris! Good catch---I had been distracted by a red herring. I appreciate your eagle eye!

Now all tests are passing, but the travis build is still exiting with 1. No errors or failures reported in log file---it successfully completes everything, then exits with 1. I've seen this one before too, and can't remember what the workaround or remedy was. Any ideas?

@ChristopherHogan @oskooi

@HomerReid
Copy link
Contributor Author

For now I reverted the .travis.yml to the version from the master repository and removed my array-slice unit test from the libmeepgeom test suite. The travis tests are now passing, verifying at least that the array-slice code successfully builds and does not screw up any existing functionality.

I have learned that by installing the 'docker' system I can recreate the travis build+test environment on my local machine, which will make it much faster to debug travis problems. I am working on setting that up. When I have made progress debugging travis builds I will add an array-slice unit test to the travis build, but it seems like that could be a separate issue that doesn't need to be a blocker for using the array-slice feature.

@oskooi @ChristopherHogan

double Scale = fmax( fabs(d1[n]), fabs(d2[n]) );
if (Scale<1.0e-8)
continue;
double RelErr = fabs(d1[n] - d2[n]) / Scale;
Copy link
Collaborator

@stevengj stevengj Sep 28, 2017

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

A more reliable way to compare two floating-point arrays is typically to compute ‖d₁-d₂‖ < reltol * (‖d₁‖+‖d₂‖)/2, where ‖⋯‖ denotes the norm (e.g. the Euclidean/L₂ norm) of the whole array.

This avoids problems with scaling and having to arbitrarily drop values below some threshold. If you want a more strict test that is sensitive to isolated errors in individual values, just use the Chebyshev (L∞) norm for ‖⋯‖.

TESTS = cyl-ellipsoid-ll

noinst_PROGRAMS = bend-flux-ll
noinst_PROGRAMS = bend-flux-ll array-slice-ll
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You're not actually running this test yet in make check, is that right?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

That's right.

// Typemap suite for array_slice
// TODO: add (cdouble *, int) version
%apply (double* INPLACE_ARRAY1, int DIM1) {(double *slice, int slice_length)};
%apply int INPLACE_ARRAY1[ANY] { int [3] };
Copy link
Collaborator

@stevengj stevengj Sep 28, 2017

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@ChristopherHogan will need to add a higher-level API function (added to the Simulation class, probably) that does the reshape etc automatically, but this is fine for now.

/***************************************************************/
#define BUFSIZE 1<<16 // use 64k buffer
if (has_imag)
{ cdouble buffer[BUFSIZE];
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

In general I would prefer not to stack-allocate here. This is allocating 2¹⁶ * (16 bytes per cdouble) = 1MB for the buffer, which could easily overflow the stack size limit on some systems. (e.g. the default stack limit on Windows is 1MB)

Use new/delete instead.

@HomerReid
Copy link
Contributor Author

HomerReid commented Sep 29, 2017

Fixed errors. Travis build/test now succeeding including libmeepgeom/array-test-ll testing array-slice functionality.
@ChristopherHogan @oskooi

@stevengj stevengj merged commit c9709ee into NanoComp:master Sep 29, 2017
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

3 participants